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Abstract 

This paper investigates the hft force exerted on an elhptical obstacle immersed in a granular flow 
through analytical calculations and computer simulations. The results are shown as a function of 
the obstacle size, orientation with respect to the flow direction (tilt angle), the restitution coefficient 
and ellipse eccentricity. The theoretical argument, based on the force exerted on the obstacle due 
to inelastic, frictionless collisions of a very dilute flow, captures the qualitative features of the lift, 
but fails to reproduce the data quantitatively. The reason behind this disagreement is that the 
dilute flow assumption on which this argument is built breaks down as a granular shock wave forms 
in front of the obstacle. More speciflcally, the shock wave change the grains impact velocity at the 
obstacle, decreasing the overall net lift obtained from a very dilute flow. 

PACS numbers: 05.10.-a, 64.70.ps, 45.70.-n, 45.70.Mg 
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I. INTRODUCTION 



Granular matter is a generic name given to a system composed of macroscopic, athermal 
particles that have mutual repulsive, dissipative interactions [1]. It is an intensely studied 
field in the Physics community given the several distinct behaviors shown by such systems 
as a consequence of different external conditions imposed on them. 

One of such conditions is that which imposes a fiow of particles, named granular fiow {2 5|. 
Within the several granular fiow examples, the fiow around immersed obstacles has received 
some attention lately 6|-[8|. One of the objectives of such investigations is to measure the 
: brce in the obstacle due to interactions with the fiowing grains, the so called granular drag 
9l-|l5|. analogously to the viscous fiow force on an obstacle. 

On the other hand, one knows that a viscous fiuid fiow around an obstacle produces an 
additional force called lift, which is perpendicular to the fiow. Given the analogy between 
a viscous fiow and the granular fiow, a lift force should exist in an obstacle immersed in 
a granular fiow under suitable conditions. However, most investigations focus only on the 
drag, while the hft studies are restricted o,dy to a few experiments and simulations QQ. 

Soller and Koehler [1^ showed that the lift on a rotating vane inserted in a granular 
packing scaled with geometric parameters of the system, such as an effective aspect ratio and 
immersion depth. The word effective means that the referred quantity, say the immersion 
depth, is considered taking into account the finite grain size. Ding et al. 15(|, dragging an 
intruder at constant velocity through a granular packing, showed that a lift force is induced 
on the intruder and it depends on its geometry. Common to both investigations is the 
fact that the lift arises as a consequence of the hydrostatic nature of the stress in granular 
packings. 

Given that these experiments were carried on such hydrostatic stress systems, variations 
of this condition might reveal different aspects of this force. For instance, is there a lift 
force in intruders immersed in flows where the stresses are not hydrostatic? If so, in what 
conditions this force arises? 

This paper is aimed at studying the lift force on an obstacle due to a dilute granular flow 
through numerical simulations. The approach will be identical to the one used in probing 



the drag on a cylinder due to a dilute granular flow 



12| . The argument drawn here to obtain 



an analytical expression for the lift shows that there is no net lift in such dilute conditions 
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on a circular obstacle. Hence, the obstacle chosen here is an elliptical one. The dependence 
of this force on the flow parameters will be obtained and studied numerically. 

Section [III is reserved for developing the theoretical argument and presenting its predic- 
tions. In section IIIIl the simulation is described along with the numerical results for the 
lift. Section HVl holds analyses regarding the argument given in the previous section in order 
to understand the discrepancies between the theoretical and the numerical results. Finally, 
section |V] has the conclusions. 

II. THEORY 

In this section, the theoretical argument leading to an expression for the lift force on 
the ellipse is developed. Also, some predictions are shown in order to be compared to the 
numerical results in the next section. 



A. Argument development 
In 



12| . the line of thought that led to an expression for the drag force on the circular 
obstacle was based on the dissipative collisions among grains and obstacle. The same line of 
thought can be drawn here in order to obtain an expression for the lift force on the ellipse 
due to collisions with the incoming stream. 

First of all, all disks are assumed to have only the horizontal velocity component, U. 
Therefore, for the ellipse given in fig. [H which is located at {xeiUe)-, a collision will occur 
only if the impact parameter (the vertical distance between a disk's center and the horizontal 
line through the ellipse's center) is in the range \yc + | ,?/d + f ] • This is so because C and 
D are, respectively, the lowest and the highest points in the ellipse, in the same way as A 
and B are the leftmost and the rightmost ones. Any collisions that happen in the segment 



AD will exert a downward lift, while collisions that occur in the AC segment will produce 



an upward lift (from now on, the AD segment, and all other segments in the text, will 



be referred to only as AD). Moreover, since the lengths of both segments are, in general. 



unequal, the longest of them, in this case AD, will suffer more collisions than the other, 
which implies a net, negative, lift force. It is this mechanism that prevents any net lift 
to take place on a circular obstacle under these conditions (in general, any body that is 
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FIG. 1: Ellipse's local frame of reference (primed system), rotated an angle 4> about its center. 
The major and minor half-axes are given by a and 6, respectively. The force F is the total (normal 
plus tangential) force between a disk and the obstacle. See text for the meaning of the four points 
A, B, C and D. 

symmetrical with respect to the flow direction, and is fully immersed in the flow, does not 
sufi^er a net lift). Therefore, the determination of the coordinates of these four points in the 
ellipse. A, C and D, is very important. 

In order to do this, one should notice that the tangent line to the ellipse is horizontal at 
D and C and vertical at A and B. Therefore, starting from the tilted eUipse's equation in 
the global frame: 



(a^ sin^ + 6^ cos^ (f) + (a^ cos^ + 6^ sin^ 0) _ 20 (a^ - h^) XY = a^6^ (1) 



where X = x — xe and Y = y — i/e, one can write for D and C: 

sin cos 0(a^ — 




(2) 



a/ sin^ (j) + b'^ cos^ 



and 




(3) 
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where the plus(minus) sign is for point D{C). A similar calculation yields, for points A and 
B: 



xb,a = \J a? cos^ (p + IP' sin 0, (4) 

and 

sin cos 0(a^ — 6^) , , 

VBA = yE± —=======, (5) 

V cos"^ + 0^ sm 

where in both eqs. the plus(minus) sign is for point B{A). 

The force due to a collision between a disk and the ellipse can be obtained by calculating 

the change in linear momentum of the disk. The total force is obtained by integrating the 

product of this individual momentum change and a suitable collision rate over the obstacle 



section facing the flow, i.e., integrating over CD. 



Suppose a disk hits the obstacle in AD. Figure |2] presents schematically the geometry of 
the collision (assumed frictionless). 




FIG. 2: Geometry of the collision between a disk and the ellipse, centered at E. They touch each 
other at the point P in the ellipse. 9 and a are the impact and collision angles. 

The grain velocity is vq = —Ui, and the impact angle is 6, which is the one between the 
direction i and the normal at the contact point. The post-collisional velocity is v = Vxi + Vyj. 
In order to obtain the final components as a function of U, one has two equations that relate 
the normal and tangential velocities before and after the collision: v ■ fi = — evo • n and 
V ■ s = Vq ■ s, where n = cos + sin 6'j, s = — sin 6i + cos are the normal and tangential 
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unit vectors at the collision point and e is the normal restitution coefficient, assumed velocity 
independent. 

From these considerations, the components of the disk velocity after the collision are 
Vx = U{s\.v?6 — ecos^O) and Vy = \U{1 + e) sin 20. Therefore, the vertical momentum 
change for a general collision is given by: 

APy{e) = mvy = ^mU{l + e) sin 29. (6) 

The number of disks that strike the ellipse in a time dt is given by the number of particles 
within the area of the parallelogram with sides Udt and ds, this last quantity being the arc 
length around the collision point. Given the geometry of the collision (figure [2]), the area of 
this parallelogram is dA = U cos Odtds. Hence: 

4z/ 4z/ 

dN = -^dA = —U cos Odtds, (7) 
ird^ ird^ 

where u is the area fraction of the incoming flow. The collision frequency is simply this 
number divided by dt. 

The arc length ds is given by the product of the collision radius, R{a), by the arc element, 
da. This product is readily evaluated for a circle. However, in the ellipse case, this is not 
so simple, since this radius varies with a. Besides, the grain size introduces an additional 
complication for calculating R{a), as seen in fig. [3], which shows the triangle formed by 
the disk center, the ellipse center and the collision point. Hence, the ellipse is assumed 
much larger than the disks and the disk size only enters the expression for R{a) as a minor 
correction. With this assumption, the arc length is given by ds = D{a)da. In parametric 
form, X = a cos a and y = 6 sin a, this length reads: 



' f dx\^ f dy^ ^ 



ds = \l j + j '^^ ~ {y cos^ a + 6^ sin a)da. 

Finally, the effective collision radius is given by: 

D{a) = (a^ cos^ a + b'^ sin^ a)^/^ (8) 

Another difficulty to obtain the expression for the lift force is that the momentum change 
(i.e., the force) depends explicitly on 9, while the collision frequency depends on a. There- 
fore, any hope of integrating the lift force over all collisions should pass through obtaining 
the relationship between the angles 6 and a. This is done in the following paragraph. 




FIG. 3: Triangle formed by the disk center, (xcyc), the elHpse center, {xe-iVe) and the colhsion 
point, {x,y) (angles are exaggerated). The angle between D{a) and the normal to the disk center 
is TT — a + 6*. Dashed lines are sketches of the surfaces of the disk and the ellipse. 

Since dx/dy gives the tangent of the angle between the normal line to a point in the 
ellipse and the x axis, from the ellipse equation ([1]), one has: 

, 2C2Y - C3X 

where Ci, C2 and C3 are the coefficients of the terms X^, Y'^ and XY in eq. ([1]). Since 
tana = Y/X and from eqs. ([2]), ([3]), (j4]) and ([5]), the relation between 6 and a can be cast 
in terms of the coordinates of the points D and A, since: 

Vd — Ve 0? sin^ + 6^ cos^ 2Ci 
" XD-XE " (l/2)sin20(a2-62) = 

and 

I/A - l/i? (1/2) sin cos 0(a2 - 6^) Cg 

tana^i 



xa — a;^ cos^ + 6^ sin 2C; 

Hence: 



1 tan — tan a .^^ 

tant^ = . (9j 

tan a A tan ao — tan a 



The lift force can now be evaluated as: 

L = - j ^PydN 

where dN = dN/dt. By using eqs. (E]), dZ]) and ([H]), the total lift on the ellipse due to 
collisions is given by: 

D 

L = -ipf/2(l + e) j sm29 cos 9ds{a), (10) 
c 

where p = Amu/nd'^ is the flow mass density and the arc length is to be calculated bearing 
in mind (jH]). 

B. Model results 

The lift force scales with the obstacle size a, as expected, since it depends linearly on 
the arc length. Also, it is seen that, for more inelastic grains the lift is smaller, since more 
inelasticity implies less intense change of momentum, see eq. (j6]). The dependence of the 
lift on the tilt angle is hidden in the integral over CD, because the lengths of AC and AD, 
which contribute forces with distinct signs, are unequal for a general 0. 

Apart from the dependence of L in these parameters, it also depends on the ellipse 
eccentricity, i.e., on the ratio of the ellipse's axes < k = b/a < 1. For k = 1, the obstacle 
is a circle, and the lift vanishes by symmetry. In the other extreme value, k = 0, which 
corresponds to a fiat plate, the lift vanishes only for = (horizontal plate) and (p = ±vr/2 
(vertical plate), which are the only symmetrical orientations with respect to an horizontal 
flow. In fig. m it is shown the predictions of eq. (fTOl) for distinct values of k (the parameter 
values in plotting these curves were chosen to agree with those used in the simulations). 

One can see that the maximum (absolute) value of L increases when k decreases. This 
happens because when A; — )■ 0, the points D and C merge with points B and A, as can be 
inferred from eqs. ([2]), ([3]), (j4]), and ([5]). In this case, AD 2a and AC 0. Since the 
lift is the difference of the contributions from AD and AC, all collisions will give positive 
contributions to the net lift value, and it should reach its maximum for a particular (p. 

The objective of the simulations is to study the lift force on the obstacle as a function 
of four parameters that appear on eq. (ITU]) , namely, the obstacle size, a, the restitution 
coefficient, e, the ellipse eccentricity, k, and the tilt angle, 0. 



8 




FIG. 4: Eq. (jlOp for distinct eccentricity, k, values. The curves were obtained with u = 0.196, 
d = 1, m = 1, [/ = 10 and e = 0.952. 

III. SIMULATION 

In this section, the simulation is detailed, along with the parameter values used in the 
computations, and the numerical results for the lift are shown. 



A. Design 

The system is composed of soft disks, with unity diameter d and unity mass m, located 
in a workspace of lengths Lx and Ly in the horizontal and vertical directions, respectively. 



They interact through normal and tangential forces, according to the model used in 13l. Il6|. 
The normal force between two disks is: 

FiJ = % + ^-, (11) 

where the first term is the conservative part, given by a simple harmonic spring force 

_ / dj + dj _ \ ^ 



where k is the spring constant, di is the i-th disk diameter, r^j is the distance between the 
two disks and fjj is a unit vector along the normal between the disks' centers. The second 
term is a dissipative, velocity dependent force, given by 

where 7^ is the normal damping coefficient and Vjj = Vj — is the relative velocity between 
the contacting disks. The tangential force at the contact point is given by the following 
expression: 

Fg = -min(7,i;,§,/iF^)vg, (12) 

where 7^ is the sliding friction constant and /i, the static friction coefficient, while v^^- is the 
unit vector along the direction of the relative velocity at the contact point. This vector is 
calculated as: 



s_ \- ( + 

where cjj is the i-th disk angular velocity and vfj = |v^ |. The total contact force Fjj = 
+ F^- vanishes if the disks are not in contact, i.e., if ^li^^ > nj. 

The values of the elastic parameters used were: k = 50000, 7^ = 10 and 7^ = 100, 
7. = 0.1 and , = 1. The restitution coefficient Q for these parameters were e = 0.952, for 
7rf = 10, and e = 0.608, for 7^ = 100. The ellipse's major and minor half-axes are a and 
b = ka, respectively. The parameter a has the values 5, 10, 20, 30 and 40, with k = 0.8. For 
the studies of the dependence on k, size a = 10 was used, with k = 0.80, 0.50 and 0.20. All 
these cases were studied for both e values. For the two smaller obstacle sizes, results were 
obtained with a. N = 2500 and Lx = Ly = 100 packing, while for the others, a N = 10000 
and Lx = Ly = 200 packing was used. In both cases, the packing fraction was about 0.200. 
The ellipse is located in xe = 0.5Lx and i/e = 0.5Ly, and its major half-axis is tilted related 
to the horizontal axis (flow direction) by (p. The values of the tilt angle were in the range 
[— 7r/2,0], divided in 7r/20 increments, which gives a total of 11 distinct values. 

At the beginning, all disks are randomly generated without overlap among them and 
the ellipse. At this stage, the obstacle is modeled as a circle with diameter 2a since this 
facilitates the overlap check. All disks have initial velocity vq = —Ui, where i is the unit 
vector in the horizontal (flow) direction and U = 10. 

The system has periodic boundaries perpendicular to the flow direction, while, along the 
flow, the conditions are the same as in 12|: whenever a disk leaves the system through 
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the right boundary, it is placed at the left one in a random vertical position. In fact, the 
code searches for a position where the incoming grain overlaps with no other disk, which 
is fairly easy, giving the low density that is used here. Its velocity is set as the incoming 
flow velocity Vq plus a small random vertical component chosen uniformly in the interval 
[—fU, fU], where / = 0.10 (results obtained with / = differ little from those shown here). 

The interactions between the disks and the ellipse are the same as those given above for 
two disks. The main problem is to resolve particle-ellipse contacts. In order to do this, all 
disks are checked for overlap with a circle of diameter 2a. If one disk overlaps with the 
circle, it is mapped in the ellipse's local frame of reference (ELFR). Then, the contact point 
coordinates, in the ELFR, are calculated using the algorithm proposed in [18] . The direction 
of the normal force is along the line joining the disk center and the contact point, while the 
tangential force is perpendicular to this direction. 

Lengths, forces and time are given in units of d, nd, and ^Jm/ respectively. Equations 
of motion are integrated with a leapfrog scheme 16|, with a time step of 0.001. Each 
simulation is performed during 10^ molecular dynamics (MD) cycles for thermalization and 
10^ MD cycles for measurements. All results are averaged over 20 and 10 independent 
runs for the = 2500 and = 10000 packings, respectively. The main interest is to 
measure the force exerted on the obstacle by the stream of grains due to the collisions 
among them. Therefore, the components of this force, drag and lift, are measured at all 
cycles after thermalization. At regular intervals, the forces are recorded as averages over 
this period. This accumulation interval is 1000 MD cycles long (results were also obtained 
with 10000 MD cycles long accumulation periods and do not differ appreciably from those 
reported here). The flow velocity and particle number density fields, v(x, y) and p{x,y), 
were measured as follows: the workspace is divided in square bins of side d. At each cycle, 
all particles are mapped in an appropriate bin and its velocity components are added to the 
respective field element. Similarly, angular profiles related to the contact angle between the 
grains and the ellipse, the collision angle a, were measured. They are the collision number 
and the velocity components. Each time there is a disk-ellipse contact, the angle formed by 
the line joining the collision point and the ellipse center with the flow direction is calculated. 
Then, the above mentioned quantities are added to the appropriate profile bins. The field 
and the profiles are presented as averages over cycles and runs. 
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FIG. 5: (Color online) Lift as a function of eccentricity, k, restitution coefficient and tilt angle, 4>. 
Squares represent k = 0.80 results, inverted triangles, k = 0.50 and left triangles, k = 0.20. Open 
and filled symbols are for e = 0.952 and 0.608. All results are for a = 10. Lines are only guides to 
the eyes. 

B. Results 

In fig. [5l the numerical results for the lift as a function of the tilt angle and the eccentricity 
are shown. 

Comparing figs. S] and 0, it is clear that eq. (ITU]) captures the qualitative features of the 
lift force, even though the correspondence becomes weak for k = 0.20 and e = 0.952. Also, 
the results agree with the prediction that the lift should be higher for less inelastic flows. 
Finally, the theoretical result overestimates the numerical ones by, roughly, a factor of 2, for 
these cases. 

In fig. [HI the numerical results for the lift force as a function of the tilt angle and 
obstacle size are shown. As expected, the lift increases with the obstacle size. However, 
the lift does not grow linearly with a. Plotting these data against the obstacle size shows 
that (L) ~ a°'''^''^\ Since the results are not conclusive, it suffices to acknowledge that more 
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FIG. 6: (Color online) Lift force as a function of the tilt angle, obstacle size, o, and restitution 
coefficient e. Distinct symbols are a = 5 (black circles), a = 10 (red squares), a = 20 (green 
diamonds), o = 30 (blue triangles) and a = 40 (orange left triangles). All are for k = 0.80. Open 
and filled symbols are results for e = 0.952 and e = 0.608, respectively. Lines are only guides to 
the eyes. 

simulations are needed in order to obtain a more reliable scaling with the obstacle size. Also, 
this figure shows that the lift force for e = 0.608 packings grows faster with a than for those 
with e = 0.952. Finally, the numerical results, as seen in fig. O also are smaller than their 
theoretical counterparts (for a = 40, the factor is about 4). This shows that the difference 
between model and theory depends only a little on the tilt angle, while the obstacle size, 
inelasticity and eccentricity are the factors that accept the most on this difference. 

The reasons behind the failure to reproduce the numerical results will be analyzed in the 
next section, where the assumption that led to eq. (JTOl) will be reviewed in detail. 
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IV. ANALYSES 



There are three basic hypotheses on which the argument for the hft force in section III Al 
was built. The first one is that particle-obstacle interactions were frictionless, despite the 
force model, eq. ( IT^ . is not. Second, the arc length expression, ds = D{a)da, does not take 
into account the fact that the grains have a finite diameter d. Finally, the collision frequency 
expression ([7]) was obtained under the hypothesis that only one particle at a time hits the 
obstacle, i.e., a disk-ellipse collision does not affect the next collision, a condition met only 
in very dilute or in ideal gas flows (which is not the case here). Therefore, its applicability 
here is clearly questionable, as inferred from the system configuration shown in fig. [71 



200 
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FIG. 7: System configuration. Parameters: a = 40, k = 0.80, e = 0.952 and 



-7r/4. 



There is a dense region that begins above point A and forms a separation boundary in 
front of the ellipse. This is the typical granular shock wave gUsI, [ii|, Q that forms around 
bodies immersed in fast granular flows. Clearly, the dilute flow assumption is not valid. 
Figure El shows the clear signature of this structure in the velocity field. 

Also, this picture shows that the two branches of the shock wave meet behind the obstacle. 
This fact could, at first sight, invalidates the lift results because this clearly does not allow 
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the use of periodic boundaries perpendicular to the stream. This is not the case, however, 
as seen in the results for the lift obtained from simulations that used rigid walls in the y 
direction. These results agree with those shown in figs. and E] within numerical error. 
Nevertheless, periodic boundaries in a particular direction should only be used when the 
system is independent in that direction. 

Before discussing the influence of the shock wave in the results, a few brief comments will 
be made regarding the first two assumptions in the theoretical argument. 

To take into account the grain size in computing the collision radius is to consider a larger 
(effective) obstacle facing the flow, which would increase the theoretical prediction for the 
lift. Since this quantity is already larger than the numerical results, it was safe to ignore it 
in the calculations. 



The influence of the friction force on the lift can be inferred from the velocity field resu 
fig. m As seen in this figure, particles follow tangential trajectories along the ellipse 



ts. 



which means that they slide along the obstacle. Since friction is a tangential force, one 



concludes that a grain sliding along AD exerts an upward lift. By symmetry, a disk exerts a 



downward lift while sliding along AC. Therefore, the net effect of friction is to decrease the 
lift exerted by collisions in each segment. The combination of both effects might decrease 
or increase the net lift. The numerical results for the lift rising from friction show that it is 
opposite to the one resulting from the normal force. In other words, friction decreases the net 
lift. The effect is small, though, because the sliding friction is also small, 7^ = 0.1. In most 
simulations of granular matter, regardless of the tangential force model, the sliding friction 
is only mildly lower than the normal interaction parameter. Here, a much lower value was 
used, which could greatly affect the results. Simulations performed with a distinct set of 
parameters, 7^ = 5000 and yU = 0.1, values which are more common to general simulations of 
granular matter, show that the results do not change significantly, and are still qualitatively 
the same as those in figs. [5]and[6l In fact, the lift data are reduced only by 5% to 10%. This 
happens due to the fact that the velocities involved in the simulations are high enough for 
the tangential force to reach the Coulomb static friction criterion, which caps its maximum 
value and limits its influence on the lift. This conclusion is supported by the fact that the 
fraction of all disk-obstacle contacts, for 7^ = 5000, that reach the failure condition is about 
99.5%, with only a very small (f) dependence. 

As stated earlier, the dilute flow assumption does not hold. A rigorous calculation of 
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the force in the obstacle, in which all dense packing effects are taken into account, is not a 



simple matter. In 



20| an attempt was made in this direction. The formation of the shock 



wave produces a dense region that shields the obstacle from the incoming particles. Since 
the interactions are dissipative, particles should reach the obstacle with reduced velocities. 
This clearly reduces the force exerted by the disks. Also, given the large density within the 
shock wave, the collision rate might be affected. Finally, as seen from the velocity field data, 
fig- El particles should hit the obstacle, on average, with a non vanishing vertical velocity 
component, i.e., particles suffer oblique inelastic collisions, which could affect the force on the 
obstacle. Both effects will be discussed in more detail in the following two subsections. This 
discussion will by made only qualitatively, since more simulations are needed to determine 
the amount of influence each of them has on the results. 



A. Obstacle shielding 

One of the consequences of the existence of the shock wave is that the obstacle is shielded 
from the flow, in a way that the grains hit it with lower horizontal velocity than the upflow 
value. This fact alone indicates that the numerical net lift should be lower than its theoretical 
prediction. Aside from reducing the incoming flow velocity, the shock wave forms a dense 
region around the obstacle through which grains should go through in order to reach it. This 
could affect the collision rate and, in turn, the net lift. Both effects must play a role in order 
to explain the faster growth of the lift force with the obstacle size for stronger inelasticity 
flow compared to that of lower inelasticity one. 

The results for the collision proflles show that, for e = 0.608, the shock wave is more 
localized around point A (where, for the k = 0.80 obstacles the peak of the proflles is 
located) and that the overall collision number is greatly increased compared to the number 
when e = 0.952. This is a consequence of the aggregation that grains suffer due to the 
inelastic interactions and is common to all cases studied. This fact implies that the force on 
the obstacle should increase due to this increase in the collision rate and could compensate, 
in part, for the decrease in the force due to the reduced incoming velocity. 

The collision proflles are changed, in a very different way, for obstacles with distinct 
eccentricities, k. Fig. IHlhas results that illustrate this effect. The most striking feature of 
these curves is that the proflles for k = 0.20 develop two peaks, instead of the more familiar 
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FIG. 9: Angular collision profiles, for a = 10, e = 0.608, k = 0.80 (upper row) and k = 0.20 (lower 
row), for three distinct tilt angles (shown in the legend box). The horizontal scales are the same 
for both plots. 

one around point A, which is still present, but it moves away from this point as the ellipse 
is oriented vertically. The other one develops around point D and is seen at tilt angles as 
high as = — tt/IO. Such features are also seen in the k = 0.50 results, although with 
smaller peaks. These results can be explained by the fact that in some cases, for stream 
velocities below a critical value, there is the ap pea rance of a gap, filled with a hot granular 
gas, between the obstacle and the shock wave HI]. In the present case, the appearance of 
the gap is a function of k alone, since the stream velocity is the same for all simulations. 
Figure [To] has a result for the density field that indicates the presence of the gap. Notice the 
denser arc formed right in front of the ellipse, the region limited by this arc and the obstacle 
is the gap. 

As a last note, the collision rate increases faster with for k = 0.20 than for k = 0.80. 
Fig. [H] has two configuration snapshots that illustrate the shock wave for two obstacles 
with distinct eccentricities. 
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FIG. 10: (Color online) Density field for a = 10, e = 0.608, k = 0.80 and (/) = -tt/2. The gap is 
the yellow (white) region between the red (light gray) arch and the ellipse (not shown). 

One can infer from these pictures the reason behind the asymmetrical peaks in the colli- 
sion profiles when = — 7r/4 and k = 0.20. There is a very small concentration of particles 
in front of AC, which certainly allows for more collisions to occur. 

These considerations still left unanswered the question about the faster growth of (L) 
with a for stronger inelasticity compared to those at lower inelasticity. Since the net lift is 
the difference of the forces exerted in AC and AD, one should seek a relative measure of 
the speeds in each of them in order to account for the size of the contributions each has 
in the net value. This is done by calculating the ratio between the average horizontal disk 
velocities in AD and AC. If this ratio is greater than 1, the collisions on AD exert, on 
average, a stronger force than those on AC. These averages are evaluated directly from the 
simulations. The results, for all tilt angles and obstacle sizes, are given in fig. [T21 

First it can be seen that except in a few cases, markedly at small obstacles, all ratios are 
larger than 1, and those for e = 0.608 cases are larger than those for 0.952. It is possible 
to identify a trend in these graphs, despite the noise: the ratio grows with obstacle size. 
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FIG. 11: Two packing configurations for a = 10, (p = — 7r/4, e = 0.608 and k = 0.80 (top) and 
k = 0.20 (bottom). 
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FIG. 12: Ratio between the average horizontal velocity components in AD and AC for obstacles 
with k = 0.80. Full lines correspond to packings with e = 0.952, while dashed lines, to e = 0.608. 
Scales in all plots are the same. 

This is consistent with the results in fig. Ej where the net lift grows faster with a for the 
e = 0.608 flows compared to those with e = 0.952. Fig. [13] has the data for the same 
quantity measured for obstacles with distinct eccentricities. 

It is seen that the ratio increases as k decreases (it reaches a factor of 2 for k = 0.20). 
These data, together with those for the collision profiles, show that the large lift observed 



for low k obstacles is a consequence not only of a larger number of collisions at AD, but 



also due to collisions with larger speeds than those at AC. 



B. Oblique impacts 

The last piece of this analysis is the role of the vertical velocity component in the lift 
value. Before proceeding, the collision argument of section [TTl is reviewed by allowing the 
incoming disks to have a vertical velocity, V, and calculating the new momentum change in 
each collision. This will elucidate the effect of the vertical velocity on the lift value. 
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FIG. 13: Ratio between the average horizontal velocity components in AD and AC for obstacles 
with k = 0.80 (left panel), k = 0.50 (middle panel) and k = 0.20 (right panel). Full lines correspond 
to packings with e = 0.952, while dashed lines, to e = 0.608. 

The velocity of a disk which collides with the obstacle is vq = —Ui + Vj. The final disk 
velocity for a inelastic, frictionless collision is obtained as: 



= U{e cos^ e - sin^ 9) -\v{l + e) sin 29 



(13) 



and 



Vy = hj{l + e) siYi29 + V{cos^9-esin^9). 



(14) 



The vertical momentum change due to this collision is given by: 

/\Py = m{vy - V), 

using the result for Vy, (fT4|) . and performing some algebra, the new momentum change is: 

"1 



APy = m(l + e) 



-U sm29 - V sin^ 9 



(15) 



where U > 0. 
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FIG. 14: Ratio between the average vertical velocity components in AD and AC. Full lines 
correspond to packings with e = 0.952, while dashed lines, to e = 0.608. Scales in all plots are the 



same. 



By looking at eq. f|T5|) one can see that a collision that happens in the range — 7r/2 < 6 < 



0, which corresponds to AC, will produce a lower momentum change if < 0. Similarly, 



for collisions that take place in the range < 9 < -k/2, which is AD, a lower momentum 
change will occur if > 0. From the velocity field data, fig. [HI it is clear that the average 



vertical velocities in AD and AC are positive and negative, respectively. The oblique impacts 



decrease the positive lift exerted by the flow in AC, but also decrease the negative lift in 



AD. Therefore, it is not obvious if both effects combined increase or decrease the net lift. 



In order for this combination of effects to increase the net lift, the vertical velocity in AD 



should be lower than the corresponding velocity in AC. In figs. [14] and [151 the absolute 
value of the ratio of the average vertical velocities in both segments is shown for distinct 
obstacle sizes eccentricities, respectively. 

It is seen that the ratio, for e = 0.608, depends little on the obstacle size and restitution 
coefficient, and that those for e = 0.952 grow with a. It varies strongly with tile angle and 
eccentricity. Despite a few cases, the ratio is larger than 1 which implies that the net lift 
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FIG. 15: Ratio between the average vertical velocity components in AD and AC. Full lines 
correspond to packings with e = 0.952, while dashed lines, to e = 0.608. The left panel is for 
k = 0.80, the middle one, for k = 0.50 and the right one, for k = 0.20. 



decreases due to oblique impacts, as provided by the existence of the shock wave. 



V. CONCLUSIONS 



A theoretical argument and numerical results on the net lift force exerted by a granular 
stream of equal disks on an ellipse were presented. The argument used to obtain a theoretical 
expression for the lift relied on inelastic, frictionless collisions of a very dilute granular flow, 
in the horizontal direction, and the ellipse. The expression obtained, eq. (fTO|) . captures 
nicely the qualitative features of the numerical results. It does not, however, reproduce the 
results quantitatively, neither does it reproduce the difference in the lift force for e = 0.952 
and e = 0.608. 

The numerical data, shown in figs. [5]and[6|, are lower than the theoretical prediction by a 
factor that depends mostly on the obstacle size (for a = 40, eq. (ITO!) is about 4 times larger 
than the measured results). 
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Additional analyses of the flow properties showed that, from the assumptions drawn in 
the text, the dilute flow one is the most problematic. The existence of the granular shock 
wave 0-8, 11, 12 1 invalidates this hypothesis. This structure affects directly the velocity 
components of the disks at contact with the obstacle as well as the collision rate. The 
horizontal velocity at contact decreases due to dissipative collisions that occur before the 
disks reach the shock wave. Also, some of the horizontal momentum is deviated by the shock 
wave around the obstacle, introducing a vertical component whose net effect is to decrease 
the lift. 

Perspectives to this work include a deeper analysis of dense packing effects on the lift. 
This would answer questions about the scaling of the lift with obstacle size, flow speed and 
density. Such analysis would also allows for a better understanding of the mechanism behind 
the appearance of the net lift, since an obvious consequence of the presence of the shock 
wave is that any collisions that occur on its edge should transmit some linear momentum 

fi, and that is when 
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22| can be good 



through the dense region until pushing the disks closets to the obstac 
the actual force is made. Force transmission models, such as those in 
starting points for this analysis. Moreover, the effect of eccentricity on the shock wave is 
something worth investigating since, as seen here, for a constant flow speed, there should be 
some eccentricity value that allows for the appearance of the gap between the shock wave 
and the obstacle. Finally, the hydrodynamical approach to granular flow could also be tested 
in such a situation. As argued in the text, ignoring the disk size does not affect the results 
for large obstacles. Hence, if the obstacle is much larger than the disks, the predictions of 
such theory should be realized in simulations, since scale separation could be achieved, at 
least approximately. 
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